Landsat collection 2, USGS

Index¶

  • Overview
  • Setup (dask, imports, query)
  • Product definition (measurements, flags)
  • Quality layer (mask)
  • Scaling and nodata
  • Visualisation
  • Appendix

Overview¶

"Since 1972, the joint NASA/ U.S. Geological Survey Landsat series of Earth Observation satellites have continuously acquired images of the Earth’s land surface, providing uninterrupted data to help land managers and policymakers make informed decisions about natural resources and the environment." https://www.usgs.gov/landsat-missions

Landsat-5, 7, 8 and 9 collection 2 products are managed by USGS. USGS make Landsat data available via number of services, including:

  • Earth Explorer - USGS View and browser
  • ESPA - USGS One-demand processing
  • AWS (STAC) - Cloud-hosted data with STAC API
  • Google Earth Engine
  • Microsoft Planetary Computer

Data source and documentation¶

Surface Reflectance, Surface Temperature and Level-1 (top of atmosphere) products for each of Landsat-5, 7, 8 and 9 are available.

EASI Asia ODC product names (Explorer): | Name | Product | Information |--|--|--| | USGS Landsat surface reflectance | landsat5_c2l2_sr | Landsat 5 Collection 2 Level-2 Surface Reflectance Product. 30m UTM based projection | | | landsat7_c2l2_sr | Landsat 7 USGS Collection 2 Surface Reflectance, processed using LEDAPS. 30m UTM based projection | | | landsat8_c2l2_sr | Landsat 8 Collection 2 Surface Reflectance, processed using LaSRC. 30m UTM based projection | | | landsat9_c2l2_sr | Landsat 9 Collection 2 Surface Reflectance, processed using LaSRC. 30m UTM based projection | | USGS Landsat surface temperature | landsat5_c2l2_st | Landsat 5 Collection 2 Level-2 UTM Surface Temperature (ST) Product | | | landsat7_c2l2_st | Landsat 7 Collection 2 Level-2 UTM Surface Temperature (ST) Product | | | landsat8_c2l2_st | Landsat 8 Collection 2 Level-2 UTM Surface Temperature (ST) Product | | | landsat9_c2l2_st | Landsat 9 Collection 2 Level-2 UTM Surface Temperature (ST) Product | | USGS Landsat Level-1 (TOA) | landsat8_c2l1 | Landsat 8 Collection 2 Level-1 (top of atmosphere) | | | landsat9_c2l1 | Landsat 9 Collection 2 Level-1 (top of atmosphere) |

EASI pipeline¶

Landsat products are read from the AWS STAC API. The data are in COG format. Two methods are shown in this notebook, each returns an essentially equivalent xarray Dataset:

  1. Read from the STAC catalog (uses odc-stac)
  2. Read from the datacube database, which has a "cached" copy of the scenes and metadata (uses odc-tools)

Notes for using the AWS STAC API:

  • Requires requester_pays = True
  • AWS source region is us-west-2 (consider egress and latency)
  • Use EASI caching-proxy settings (to help reduce egress and latency costs)

Setup¶

In [1]:
%matplotlib inline

# Data tools
import numpy as np
import xarray as xr
import pandas as pd
from datetime import datetime
import matplotlib.pyplot as plt

# Formatting options
pd.set_option("display.max_rows", None)
# plt.rcParams['figure.figsize'] = [12, 8]

# Datacube
import datacube
from datacube.utils import masking  # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/masking.py
from odc.algo import enum_to_bool   # https://github.com/opendatacube/odc-tools/blob/develop/libs/algo/odc/algo/_masking.py
from datacube.utils.aws import configure_s3_access

# Notebook helper tools (in dea_tools or in this repo)
import sys
from os import environ
repo = f'{environ["HOME"]}/eocsi-hackathon-2022'  # No easy way to get repo directory
if repo not in sys.path: sys.path.append(repo)
from tools.notebook_utils import xarray_object_size
try:
    from dea_tools.plotting import display_map, rgb
    from dea_tools.datahandling import mostcommon_crs
except ImportError:
    # Local copy of selected dea_tools
    from tools.datacube_utils import display_map, mostcommon_crs
    rgb = None  # Not copied or adapted yet

# Holoviews, Datashader and Bokeh
import hvplot.pandas
import hvplot.xarray
import holoviews as hv
import panel as pn
import colorcet as cc
import cartopy.crs as ccrs
from datashader import reductions
from holoviews import opts
# import geoviews as gv
# from holoviews.operation.datashader import rasterize
hv.extension('bokeh', logo=False);

Connect to the ODC database¶

In [2]:
# Initialise a datacube connection
dc = datacube.Datacube()

# Access AWS "requester-pays" buckets
configure_s3_access(aws_unsigned=False, requester_pays=True)

# Use EASI caching-proxy (applies to selected source buckets)
environ["AWS_HTTPS"] = "NO"
environ["GDAL_HTTP_PROXY"] = "easi-caching-proxy.caching-proxy:80"
print(f'Will use caching proxy at: {environ.get("GDAL_HTTP_PROXY")}')
Will use caching proxy at: easi-caching-proxy.caching-proxy:80

Example query¶

In [3]:
# Vietnam - Ha Long
latitude = (20.7, 21.1)
longitude = (106.7, 107.2)
time=('2020-02-01', '2020-04-01')
display_map(longitude, latitude)
Out[3]:
Make this Notebook Trusted to load map: File -> Trust Notebook

1. Read from the STAC catalog¶

In [ ]:
 
In [ ]:
 

2. Read from the ODC database¶

In [8]:
# Select a product
product = 'landsat8_c2l2_sr'

# Split the query to determine the most common CRS (essentially call find_datasets())
query = {
    'x': longitude,    # "x" axis bounds
    'y': latitude,      # "y" axis bounds
    'time': time,           # Any parsable date strings
}

# Most common CRS
# dea_tools: ValueError if query does not return datasets
native_crs = mostcommon_crs(dc, product, query)
print(f'Native CRS: {native_crs}')

query.update({
    'product': product,                     # Product name
    'output_crs': native_crs,               # EPSG code
    'resolution': (30, 30),                # Target resolution
    'group_by': 'solar_day',                # Scene ordering
    # 'dask_chunks': {'x': 2048, 'y': 2048},  # Dask chunks
})

# Load data
data = dc.load(**query)
display(xarray_object_size(data))
display(data)
Native CRS: EPSG:32648
'Dataset size: 47.61 MB'
<xarray.Dataset>
Dimensions:      (time: 1, y: 1498, x: 1753)
Coordinates:
  * time         (time) datetime64[ns] 2020-02-14T03:17:16.453607
  * y            (y) float64 2.29e+06 2.29e+06 2.29e+06 ... 2.335e+06 2.335e+06
  * x            (x) float64 6.766e+05 6.766e+05 ... 7.291e+05 7.291e+05
    spatial_ref  int32 32648
Data variables:
    coastal      (time, y, x) uint16 0 0 0 0 0 ... 43982 44124 44485 44905 45366
    blue         (time, y, x) uint16 0 0 0 0 0 ... 43958 44162 44487 44933 45458
    green        (time, y, x) uint16 0 0 0 0 0 ... 42835 42980 43240 43892 44440
    red          (time, y, x) uint16 0 0 0 0 0 ... 42556 42741 43012 43660 44289
    nir08        (time, y, x) uint16 0 0 0 0 0 ... 41679 41879 42211 42768 43458
    swir16       (time, y, x) uint16 0 0 0 0 0 ... 29493 29956 30389 32847 33978
    swir22       (time, y, x) uint16 0 0 0 0 0 ... 24751 25260 25750 28564 29340
    qa_pixel     (time, y, x) uint16 1 1 1 1 1 ... 22280 22280 22280 22280 22280
    qa_aerosol   (time, y, x) uint8 1 1 1 1 1 1 1 ... 224 192 224 224 192 224
    qa_radsat    (time, y, x) uint16 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
Attributes:
    crs:           EPSG:32648
    grid_mapping:  spatial_ref
xarray.Dataset
    • time: 1
    • y: 1498
    • x: 1753
    • time
      (time)
      datetime64[ns]
      2020-02-14T03:17:16.453607
      units :
      seconds since 1970-01-01 00:00:00
      array(['2020-02-14T03:17:16.453607000'], dtype='datetime64[ns]')
    • y
      (y)
      float64
      2.29e+06 2.29e+06 ... 2.335e+06
      units :
      metre
      resolution :
      30.0
      crs :
      EPSG:32648
      array([2289885., 2289915., 2289945., ..., 2334735., 2334765., 2334795.])
    • x
      (x)
      float64
      6.766e+05 6.766e+05 ... 7.291e+05
      units :
      metre
      resolution :
      30.0
      crs :
      EPSG:32648
      array([676575., 676605., 676635., ..., 729075., 729105., 729135.])
    • spatial_ref
      ()
      int32
      32648
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 48N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32648"]]
      grid_mapping_name :
      transverse_mercator
      array(32648, dtype=int32)
    • coastal
      (time, y, x)
      uint16
      0 0 0 0 ... 44124 44485 44905 45366
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32648
      grid_mapping :
      spatial_ref
      array([[[    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              ...,
              [34131, 33715, 33247, ..., 45812, 46125, 46120],
              [33359, 33044, 32774, ..., 45391, 45910, 45997],
              [33305, 33116, 32975, ..., 44485, 44905, 45366]]], dtype=uint16)
    • blue
      (time, y, x)
      uint16
      0 0 0 0 ... 44162 44487 44933 45458
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32648
      grid_mapping :
      spatial_ref
      array([[[    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              ...,
              [33767, 33227, 32851, ..., 46046, 46331, 46427],
              [33095, 33026, 32690, ..., 45391, 45999, 46122],
              [33182, 32924, 32842, ..., 44487, 44933, 45458]]], dtype=uint16)
    • green
      (time, y, x)
      uint16
      0 0 0 0 ... 42980 43240 43892 44440
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32648
      grid_mapping :
      spatial_ref
      array([[[    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              ...,
              [32740, 32424, 32305, ..., 44688, 44949, 44959],
              [31851, 31409, 31029, ..., 44300, 44867, 45124],
              [31543, 31357, 31283, ..., 43240, 43892, 44440]]], dtype=uint16)
    • red
      (time, y, x)
      uint16
      0 0 0 0 ... 42741 43012 43660 44289
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32648
      grid_mapping :
      spatial_ref
      array([[[    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              ...,
              [32214, 31657, 31352, ..., 44777, 45083, 45307],
              [31166, 30948, 30664, ..., 44138, 44785, 45089],
              [30994, 30815, 30758, ..., 43012, 43660, 44289]]], dtype=uint16)
    • nir08
      (time, y, x)
      uint16
      0 0 0 0 ... 41879 42211 42768 43458
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32648
      grid_mapping :
      spatial_ref
      array([[[    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              ...,
              [31658, 31147, 30910, ..., 44213, 44509, 44818],
              [30751, 30511, 30344, ..., 43096, 43966, 44376],
              [30694, 30605, 30571, ..., 42211, 42768, 43458]]], dtype=uint16)
    • swir16
      (time, y, x)
      uint16
      0 0 0 0 ... 29956 30389 32847 33978
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32648
      grid_mapping :
      spatial_ref
      array([[[    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              ...,
              [15358, 15402, 14928, ..., 33580, 34443, 34444],
              [14230, 13606, 13296, ..., 33782, 34949, 35143],
              [13250, 12678, 12519, ..., 30389, 32847, 33978]]], dtype=uint16)
    • swir22
      (time, y, x)
      uint16
      0 0 0 0 ... 25260 25750 28564 29340
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32648
      grid_mapping :
      spatial_ref
      array([[[    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              [    0,     0,     0, ...,     0,     0,     0],
              ...,
              [11482, 11630, 11297, ..., 29348, 29948, 29727],
              [10514, 10096,  9807, ..., 29307, 30573, 30570],
              [10040,  9741,  9643, ..., 25750, 28564, 29340]]], dtype=uint16)
    • qa_pixel
      (time, y, x)
      uint16
      1 1 1 1 ... 22280 22280 22280 22280
      units :
      bit_index
      nodata :
      1
      flags_definition :
      {'snow': {'bits': 5, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'clear': {'bits': 6, 'values': {'0': 'not_clear', '1': 'clear'}}, 'cloud': {'bits': 3, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'water': {'bits': 7, 'values': {'0': 'land_or_cloud', '1': 'water'}}, 'cirrus': {'bits': 2, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'nodata': {'bits': 0, 'values': {'0': False, '1': True}}, 'qa_pixel': {'bits': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], 'values': {'1': 'Fill', '2': 'Dilated Cloud', '4': 'Cirrus', '8': 'Cloud', '16': 'Cloud Shadow', '32': 'Snow', '64': 'Clear', '128': 'Water', '256': 'Cloud Confidence low bit', '512': 'Cloud Confidence high bit', '1024': 'Cloud Shadow Confidence low bit', '2048': 'Cloud Shadow Confidence high bit', '4096': 'Snow Ice Confidence low bit', '8192': 'Snow Ice Confidence high bit', '16384': 'Cirrus Confidence low bit', '32768': 'Cirrus Confidence high bit'}, 'description': 'Level 2 pixel quality'}, 'cloud_shadow': {'bits': 4, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'dilated_cloud': {'bits': 1, 'values': {'0': 'not_dilated', '1': 'dilated'}}, 'cloud_confidence': {'bits': [8, 9], 'values': {'0': 'none', '1': 'low', '2': 'medium', '3': 'high'}}, 'cirrus_confidence': {'bits': [14, 15], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'snow_ice_confidence': {'bits': [12, 13], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'cloud_shadow_confidence': {'bits': [10, 11], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}}
      crs :
      EPSG:32648
      grid_mapping :
      spatial_ref
      array([[[    1,     1,     1, ...,     1,     1,     1],
              [    1,     1,     1, ...,     1,     1,     1],
              [    1,     1,     1, ...,     1,     1,     1],
              ...,
              [22280, 22280, 22280, ..., 22280, 22280, 22280],
              [22280, 22280, 22280, ..., 22280, 22280, 22280],
              [22280, 22280, 24082, ..., 22280, 22280, 22280]]], dtype=uint16)
    • qa_aerosol
      (time, y, x)
      uint8
      1 1 1 1 1 1 ... 192 224 224 192 224
      units :
      bit_index
      nodata :
      1
      flags_definition :
      {'water': {'bits': 2, 'values': {'0': 'not_water', '1': 'water'}}, 'nodata': {'bits': 0, 'values': {'0': False, '1': True}}, 'qa_aerosol': {'bits': [0, 1, 2, 3, 4, 5, 6, 7], 'values': {'1': 'Fill', '2': 'Valid aerosol retrieval', '4': 'Water', '8': 'Unused', '16': 'Unused', '32': 'Interpolated Aerosol', '64': 'Aerosol Level low bit', '128': 'Aerosol Level high bit'}, 'description': 'Aerosol quality assessment'}, 'aerosol_level': {'bits': [6, 7], 'values': {'0': 'climatology', '1': 'low', '2': 'medium', '3': 'high'}}, 'valid_retrieval': {'bits': 1, 'values': {'0': 'not_valid', '1': 'valid'}}, 'interp_retrieval': {'bits': 5, 'values': {'0': 'not_aerosol_interpolated', '1': 'aerosol_interpolated'}}}
      crs :
      EPSG:32648
      grid_mapping :
      spatial_ref
      array([[[  1,   1,   1, ...,   1,   1,   1],
              [  1,   1,   1, ...,   1,   1,   1],
              [  1,   1,   1, ...,   1,   1,   1],
              ...,
              [224, 224, 224, ..., 224, 224, 224],
              [224, 224, 224, ..., 224, 224, 224],
              [224, 224, 192, ..., 224, 192, 224]]], dtype=uint8)
    • qa_radsat
      (time, y, x)
      uint16
      0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0
      units :
      bit_index
      nodata :
      1
      flags_definition :
      {'qa_radsat': {'bits': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], 'values': {'1': 'Band 1 Data Saturation', '2': 'Band 2 Data Saturation', '4': 'Band 3 Data Saturation', '8': 'Band 4 Data Saturation', '16': 'Band 5 Data Saturation', '32': 'Band 6 Data Saturation', '64': 'Band 7 Data Saturation', '128': 'Unused', '256': 'Band 9 Data Saturation', '512': 'Unused', '1024': 'Unused', '2048': 'Terrain occlusion'}, 'description': 'Radiometric saturation'}, 'b1_saturation': {'bits': 0, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b2_saturation': {'bits': 1, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b3_saturation': {'bits': 2, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b4_saturation': {'bits': 3, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b5_saturation': {'bits': 4, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b6_saturation': {'bits': 5, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b7_saturation': {'bits': 6, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'b9_saturation': {'bits': 8, 'values': {'0': 'no_saturation', '1': 'saturated_data'}}, 'terrain_occlusion': {'bits': 11, 'values': {'0': 'no_terrain_occlusion', '1': 'terrain_occlusion'}}}
      crs :
      EPSG:32648
      grid_mapping :
      spatial_ref
      array([[[0, 0, 0, ..., 0, 0, 0],
              [0, 0, 0, ..., 0, 0, 0],
              [0, 0, 0, ..., 0, 0, 0],
              ...,
              [0, 0, 0, ..., 0, 0, 0],
              [0, 0, 0, ..., 0, 0, 0],
              [0, 0, 0, ..., 0, 0, 0]]], dtype=uint16)
  • crs :
    EPSG:32648
    grid_mapping :
    spatial_ref
In [ ]:
# Optional. Filter Datasets prior to Load

# For example, to load only Tier 1 (best quality) datasets and exclude Tier 2 (good quality).
# See https://www.usgs.gov/core-science-systems/nli/landsat/landsat-collection-1 for a description of Landsat processing Tiers.

# 1. Remove load parameters from query object
# Borrowed from https://github.com/GeoscienceAustralia/dea-notebooks/blob/develop/Tools/dea_tools/datahandling.py

# non_load_query = datacube_utils.dc_query_only(**query)
# dataset_list = dc.find_datasets(**non_load_query)

# 2. Check your query has results
# Borrowed from https://github.com/GeoscienceAustralia/dea-notebooks/blob/develop/Tools/dea_tools/datahandling.py

# if len(dataset_list) == 0:
#    print("No data available for query: ensure that "
#          "the products specified have data for the "
#          "time and location requested")

# 3. Check what details are available in each dataset
# See https://github.com/opendatacube/datacube-core/blob/develop/datacube/model/__init__.py, class Dataset
# display(datasets[0].__dict__)

# 4. Filter based on property of interest
# For Landsat, the Tier label is available in the 'landsat:collection_category' property

# dataset_list = [ds for ds in dataset_list if ds.metadata_doc['properties']['landsat:collection_category'] == 'T1']
# if len(dataset_list) == 0:
#    print("No data available after filtering")

# 5. Update query object for next cell
# 'datasets' will used instead of the standard database lookup

# query['datasets'] = dataset_list
# query
In [ ]:
# Load data
data = dc.load(**query)

notebook_utils.heading(notebook_utils.xarray_object_size(data))
display(data)

# Calculate valid (not nodata) masks for each layer
valid_mask = masking.valid_data_mask(data)
notebook_utils.heading('Valid data masks for each variable')
display(valid_mask)

Product definition¶

Display the measurement definitions for the selected product.

Use list_measurements to show the details for a product, and masking.describe_variable_flags to show the flag definitions.

In [ ]:
# Measurement definitions for the selected product
measurement_info = dc.list_measurements().loc[query['product']]
notebook_utils.heading(f'Measurement table for product: {query["product"]}')
notebook_utils.display_table(measurement_info)  # Default pandas table display. Some rows or columns may be abbreviated

# Separate lists of measurement names and flag names
measurement_names = measurement_info[ pd.isnull(measurement_info.flags_definition)].index
flag_names        = measurement_info[pd.notnull(measurement_info.flags_definition)].index

notebook_utils.heading('Selected Measurement and Flag names')
notebook_utils.display_table(pd.DataFrame({
    'group': ['Measurement names', 'Flag names'],
    'names': [', '.join(measurement_names), ', '.join(flag_names)]
}))

# Flag definitions
for flag in flag_names:
    notebook_utils.heading(f'Flag definition table for flag name: {flag}')
    notebook_utils.display_table(masking.describe_variable_flags(data[flag]))

Quality layer¶

In [ ]:
# Make L2_FLAGS image
flag_name = 'pixel_qa'
flag_data = data[[flag_name]].where(valid_mask[flag_name]).persist()   # Dataset
display(flag_data)
In [ ]:
##
## Drat!
## Geoviews/Cartopy projection from epsg:32649 to PlateCarree doesn't work (works for S2 UTM)
##

# These options manipulate the color map and colorbar to show the categories for this product
options = {
    'title': f'Flag data for: {query["product"]} ({flag_name})',
    'cmap': cc.rainbow,
    'colorbar': True,
    'width': 700,
    'height': 450,
    'aspect': 'equal',
    'tools': ['hover'],
}

# Set the Dataset CRS
plot_crs = native_crs
if plot_crs == 'epsg:4326':
    plot_crs = ccrs.PlateCarree()


# Native data and coastline overlay:
# - Comment `crs`, `projection`, `coastline` to plot in native_crs coords
# TODO: Update the axis labels to 'longitude', 'latitude' if `coastline` is used

quality_plot = flag_data.hvplot.image(
    x = 'x', y = 'y',         # Dataset x,y dimension names
    rasterize = True,                        # Use Datashader
    aggregator = reductions.mode(),          # Datashader selects mode value, requires 'hv.Image'
    precompute = True,                       # Datashader precomputes what it can
    crs = plot_crs,                          # Dataset CRS
    projection = ccrs.PlateCarree(),         # Output Projection (ccrs.PlateCarree() when coastline=True)
    coastline = '10m',                       # Coastline = '10m'/'50m'/'110m'
).options(opts.Image(**options)).hist()

# display(quality_plot)
# Optional: Change the default time slider to a dropdown list, https://stackoverflow.com/a/54912917
fig = pn.panel(quality_plot, widgets={'time': pn.widgets.Select})  # widget_location='top_left'
display(fig)
In [ ]:
# Create mask layer

# "L3 Mask Default"
good_pixel_flags = {
    'snow': 'no_snow',                    # 'no_snow', 'snow'
#     'clear': 'clear_land',                # 'no_clear_land', 'clear_land'
    'cloud': 'no_cloud',                  # 'no_cloud', 'cloud'
#     'water': 'water',                     # 'no_water', 'water'
    'nodata': False,                      # False, True
    'cloud_shadow': 'no_cloud_shadow',    # 'no_cloud_shadow', 'cloud_shadow'
#     'cloud_confidence': 'medium',         # 'none', 'low', 'medium', 'high'
#     'cirrus_confidence': 'medium',        # 'none', 'low', 'medium', 'high'
#     'terrain_occlusion': 'no_occlusion',  # 'no_occlusion', 'occlusion'
}

good_pixel_mask = masking.make_mask(data[flag_name], **good_pixel_flags)  # -> bool array
display(good_pixel_mask)  # Type: bool

Scaling and nodata¶

In [ ]:
# Define scaling

# usgs_espa_ls8c1_sr, usgs_espa_ls8c1_ar
scale = 0.0001
offset = 0.
In [ ]:
# Select a layer and apply masking and scaling, then persist in dask

layer_name = 'nir'

# Apply valid mask and good pixel mask
layer = data[[layer_name]].where(valid_mask[layer_name] & good_pixel_mask) * scale + offset
layer = layer.persist()

Visualisation¶

In [ ]:
##
## Drat!
## Geoviews/Cartopy projection from epsg:32649 to PlateCarree doesn't work (works for S2 UTM)
##

# Generate a plot

options = {
    'title': f'{query["product"]}: {layer_name}',
    'width': 1000,
    'height': 450,
    'aspect': 'equal',
    'cmap': cc.rainbow,
    'clim': (0, 0.05),                          # Limit the color range depending on the layer_name
    'colorbar': True,
    'tools': ['hover'],
}

# Set the Dataset CRS
plot_crs = native_crs
if plot_crs == 'epsg:4326':
    plot_crs = ccrs.PlateCarree()


# Native data and coastline overlay:
# - Comment `crs`, `projection`, `coastline` to plot in native_crs coords
# TODO: Update the axis labels to 'longitude', 'latitude' if `coastline` is used
    
layer_plot = layer.hvplot.image(
    x = 'x', y = 'y',                        # Dataset x,y dimension names
    rasterize = True,                        # Use Datashader
    aggregator = reductions.mean(),          # Datashader selects mean value
    precompute = True,                       # Datashader precomputes what it can
    crs = plot_crs,                          # Dataset crs
    projection = ccrs.PlateCarree(),         # Output projection (use ccrs.PlateCarree() when coastline=True)
    coastline='10m',                         # Coastline = '10m'/'50m'/'110m'
).options(opts.Image(**options)).hist(bin_range = options['clim'])

# display(layer_plot)
# Optional: Change the default time slider to a dropdown list, https://stackoverflow.com/a/54912917
fig = pn.panel(layer_plot, widgets={'time': pn.widgets.Select})  # widget_location='top_left'
display(fig)

# Good image of Sarawak: 2020-01-27 0555

Appendix¶

Reference material

In [ ]:
# Filter datasets

# From load_ard()
datasets = dc.find_datasets(product=product, **query)

# Remove datasets after the Landsat 7 SLC failure, May 31 2003.
if product in ('ga_ls7e_ard_3', 'usgs_espa_ls7c1_sr'):
    datasets = [i for i in datasets if
                normalise_dt(i.time.begin) <
                datetime.datetime(2003, 5, 31)]
In [ ]:
# Masking and Scaling

# dea-notebooks/Real_world_examples/ARD_Intercomparison/utilities/util.py

ls8_USGS_cloud_pixel_qa_value = [324, 352, 368, 386, 388, 392, 400, 416, 
                                     432, 480, 864, 880, 898, 900, 904, 928, 
                                     944, 992, 1350]
non_ls8_USGS_cloud_pixel_qa_value = [72, 96, 112, 130, 132, 136, 144, 160, 
                                         176, 224]
    non_ls8_USGS_sr_cloud_qa_value = [2, 4, 12, 20, 34, 36, 52]

mask_data = data[mask_band]
    nodata_value = mask_data.nodata
    nodata_cloud_value = []
    
    if 'usgs' in source_prod:
        if 'ls8' in source_prod:
            nodata_cloud_value = ls8_USGS_cloud_pixel_qa_value
        else:
            if mask_band == 'sr_cloud_qa':
                nodata_cloud_value = non_ls8_USGS_sr_cloud_qa_value
            else:
                nodata_cloud_value = non_ls8_USGS_cloud_pixel_qa_value
                
        nodata_cloud_value.append(nodata_value)
        nodata_cloud = np.isin(mask_data, nodata_cloud_value) 
        cld_free = data.where(~nodata_cloud).dropna(dim='time', how='all')
        
    # remove nodata for the pixel of interest
    cld_free_valid = masking.mask_invalid_data(cld_free)
In [ ]:
# Visualisation